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ABSTRACT 


This  report  introduces  new  research  on  the  generalised  Marcum  Q-Function, 
and  in  particular,  the  probability  of  detection  of  a  number  of  incoherently  inte¬ 
grated  signals  in  a  Gaussian  clutter  and  noise  environment.  A  new  probabilistic 
association  is  derived,  linking  this  detection  probability  with  a  probability  as¬ 
sociated  with  two  independent  Poisson  random  variables.  Additionally,  it  is 
shown  that  this  detection  probability  is  the  solution  to  two  stochastic  Yolterra 
integral  equations.  This  results  in  a  means  of  obtaining  estimates  of  this  detec¬ 
tion  probability.  Specifically,  lower  and  upper  bounds  are  derived  using  these 
representations,  and  the  bounds  are  compared  with  known  results.  As  a  by¬ 
product  of  this  work  a  new  useful  expression  for  the  differences  in  distributions 
of  independent  Poissons  random  variables  is  obtained. 
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Stochastic  Representations  of  the  Marcum  Q-Function  and 
Associated  Radar  Detection  Probabilities 


EXECUTIVE  SUMMARY 


The  work  presented  here  is  a  consequence  of  the  ongoing  long  range  research  associated 
with  radar  detection  issues  that  arose  out  of  the  task  AIR  01/217,  which  has  now  been 
succeeded  by  AIR  04/216.  The  purpose  of  this  task  is  to  provide  the  Royal  Australian  Air 
Force  with  technical  advice  on  the  performance  of  the  Elta  EL/M-2022  maritime  radar. 
The  latter  radar  is  used  in  the  AP-3C  Orion  fleet.  Key  performance  measures  of  a  radar 
include  probabilities  of  false  alarm  and  detection.  Earlier  work  by  the  author  focused 
on  how  Monte  Carlo  methods  could  be  used  to  estimate  such  quantities.  The  research 
found  in  this  report  arose  out  of  the  need  to  find  an  efficient  Monte  Carlo  estimator  for 
a  single  pulse  probability  of  detection  of  a  target  in  Gaussian  clutter  and  noise.  This 
was  a  component  in  a  two-tiered  Monte  Carlo  estimator  for  a  binary  integrated  detection 
probability.  A  new  association  was  discovered  between  the  pulse  detection  probability  and 
a  pair  of  independent  Poisson  random  variables.  The  mathematics  behind  this  discovery 
is  presented  here.  Furthermore,  it  turns  out  that  this  result  can  be  extended  to  the  case 
of  a  series  of  incoherently  integrated  pulses.  Consequently,  this  provides  a  new  twist  on 
the  generalised  Marcum  Q-Function.  The  latter  is  an  important  function  in  the  study 
of  radar  and  communications,  and  hence  new  representations  for  it  are  of  interest  to  the 
wider  radar  community. 

In  addition  to  presenting  this  new  result,  it  is  shown  that  this  detection  probability  is  the 
solution  to  two  stochastic  Volterra  integral  equations  of  the  second  kind.  These  represen¬ 
tations  for  the  detection  probability  suggest  ways  in  which  bounds  can  be  obtained.  The 
importance  of  bounds  on  such  probabilities  is  that  they  indicate  the  minimum  and  maxi¬ 
mum  performance  levels  of  a  radar  detection  scheme.  We  construct  new  lower  and  upper 
bounds  on  a  specific  detection  probability,  and  compare  them  with  well-known  results 
from  the  signal  processing  literature. 
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8  The  same  plot  as  for  Figure  7,  except  with  the  removal  of  the  bound  (59)  due 
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Glossary 

IN  Natural  numbers  {0, 1,2,.. 

P  Probability. 


E  Statistical  expectation. 


I  Indicator  function: 


1  if  x  £  A; 

0  otherwise. 


:=  Defined  to  be. 


a  V  b  Maximum  of  a  and  b. 

=  Equality  in  distribution:  X  =  Y  is  equivalent  to  P(X  e4)  =  P(F  €  A)  for  all  sets  A. 
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Po(A)  Poisson  Distribution  with  mean  A  >  0:  if  X  =  Po(A),  then  P(X  =  j)  =  e  j*J ,  for 
all  j  €  IN. 

Po(A){^4}  Cumulative  Poisson  probability  on  set  iclN:  Po(A){Al}  =  J2jeA  e  jr~  ■ 

Bin(n,p)  Binomial  Distribution  with  parameters  n  e  IN  and  p  €  [0, 1]:  if  X  =  Bin(n,p), 
then  P(X  =  j)  =  (1  —  p)n~^ ,  for  j  €  {0, 1, . . . ,  n}. 

R(a,/3)  Uniform  (or  rectangular)  distribution  on  the  interval  [a,/3\  (a  <  (3):  If  X  = 
R(a,/3),  then  P(X  <  x)  =  f°r  x  €  [ct,/3]. 
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1  Preliminaries 


1.1  Background 


The  single  pulse  probability  of  detection,  for  a  constant  target  in  Gaussian  clutter  and 
noise,  has  been  the  subject  of  much  investigation  in  foundational  studies  in  the  mathemat¬ 
ical  analysis  of  radar  detectors.  Expressions  for  such  probabilities  first  appeared  in  [Rice 
1944  and  1945,  Marcum  1960  and  Di  Franco  and  Rubin  1968],  and  more  recently,  [Lev- 
anon  1988  and  Minkler  and  Minkler  1990]  contain  detailed  analyses  deriving  this  classical 
detection  probability. 

Under  a  Neyman-Pearson  regime,  this  probability  of  detection  appears  as  an  integral  in¬ 
volving  a  modified  Bessel  function  of  order  zero.  This  integral  does  not  have  a  closed 
analytic  form.  As  it  appears  in  a  number  of  modelling  applications,  such  as  binary  in¬ 
tegration  [Shnidman  1998],  a  number  of  authors  have  examined  ways  of  estimating  it 
efficiently.  [Shnidman  1995,  1998  and  2002]  applied  a  Maclaurin  series  expansion  to  the 
Bessel  function,  and  consequently  obtains  a  number  of  useful  results.  Simulation  method¬ 
ology  can  also  be  used  to  estimate  this  pulse  probability  of  detection.  In  [Weinberg  and 
Kyprianou  2005],  Monte  Carlo  methods  are  used  to  estimate  it. 

As  a  result  of  searching  for  an  efficient  estimator  for  this  pulse  detection  probability,  a  new 
probabilistic  expression  for  it  has  been  obtained.  It  can  be  shown  [Weinberg  and  Kyprianou 
2005]  that  the  detection  probability  is  identical  to  a  comparison  of  two  independent  Poisson 
random  variables.  One  is  centred  on  the  average  signal  strength,  while  the  other  is  centred 
on  a  normalised  detection  threshold. 


1.2  Probability  of  Detection  in  Gaussian  Clutter 


The  following  is  included  for  completeness,  and  to  provide  a  context  for  the  work  to  be 
discussed  in  subsequent  sections,  and  is  based  upon  the  developments  in  [Levanon  1988]. 
Consider  a  radar  operating  in  Gaussian  clutter  and  noise,  from  which  a  single  pulse  is 
transmitted.  The  transmitted  signal  is  a  sine  wave  with  period  if)  and  frequency  oj.  The 
returned  signal  will  be  assumed  to  be  a  phase  shifted  version  of  the  original,  with  the 
addition  of  interference.  For  modelling  simplicity,  we  do  not  differentiate  between  clutter 
and  noise,  nor  any  other  environmental  factors  that  may  distort  the  signal.  We  will 
just  assume  the  total  interference  is  a  Gaussian  random  variable.  The  returned  signal  is 
passed  to  a  narrow  bandpass  filter,  with  centre  frequency  u.  We  assume  this  filter  has 
a  rectangular  response  with  bandwidth  fs-  Then  assuming  that  /b  >  ^,  the  returned 
signal  is 

s(t)  =  Acos{ut  —  6)  =  acos(W)  +  bsm(u>t),  (1) 

where  0  =  arctan  is  the  phase  shift  of  the  signal,  and  the  amplitude  is  A  =  Va2  +  b2. 
We  assume  that  a  and  b,  and  consequently  A,  are  deterministic,  and  that  6  is  uniformly 
distributed  on  the  interval  [0,27 r). 
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When  Gaussian  noise  is  passed  through  a  narrow  bandpass  filter,  the  output  can  be  written 
as 

n(t)  =  X(t)  cos(c jjt)  +  Y ( t )  sin(u>t),  (2) 

where  X(t)  and  Y(t)  are  both  independent  and  identically  distributed  Gaussian  random 
variables,  with  zero  mean  and  variance  a2.  By  combining  both  (1)  and  (2),  the  radar 
signal  return,  in  additive  Gaussian  noise,  at  the  detector  is 

C(t)  =  s(t)  +  n(t ) 

=  (A  cos#  +  X(t))  cos(cuf)  +  (A  sin#  +  Y(t))  sin(aA) 

=  (a  +  X (t))  cos(ut)  +  (b  +  Y (t))  sm(u}t)  (3) 


where  Iq,  the  modified  Bessel  function,  of  the  first  kind,  of  order  zero,  is  defined  to  be 

1  /*7r 

h{x)=  exp(xcos  (9-ip))d9=—  excos6de.  (7) 

Jo  27T  J—tt 

Note  that  (6)  does  not  depend  on  9,  so  that  fii(r\0)  =  /i?(r).  We  are  now  in  a  position 
where  we  can  specify  the  probabilities  of  false  alarm  and  detection  for  a  single  pulse.  The 
classical  approach  to  radar  detection  is  to  declare  a  target  present  when  some  statistic 
of  the  returned  signal  exceeds  some  pre-dehned  threshold  level.  A  common  choice  in  the 
literature  is  for  the  amplitude  ( R(t ))  to  exceed  a  mean  level  [see  Levanon  1988].  We 
suppose  this  threshold  is  nT.  Observe  that,  in  view  of  (6),  in  the  case  where  there  is  no 
signal  present  in  the  return,  so  that  A  =  0,  the  pdf  of  the  amplitude  is  Rayleigh  distributed 
with  parameter  a.  The  probability  of  false  alarm  is  the  probability  of  declaring  a  target 
present  when  there  is  only  noise.  Hence,  under  the  Neyman-Pearson  criterion, 
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The  probability  of  detection  is  the  probability  of  correctly  declaring  a  target  present.  Un¬ 
der  the  Neyman-Pearson  Criterion,  for  the  case  of  a  band-limited  signal  (1)  in  bandlimited 
Gaussian  noise,  it  is  given  by  the  integral 


POO 

r 2  +  A2 

To 

rA' 

/  o 

2a2 

.  v2 . 

dr, 


(9) 


where  the  SNR  is  given  by  =  ^4  [see  Levanon  1988].  The  larger  the  signal  to  noise 
strength,  the  more  likely  we  are  to  detect  the  target.  In  order  to  emphasize  the  dependence 
of  (9)  on  the  SNR  q  and  the  threshold  r,  we  will  write  the  detection  probability  as  p(s,  r). 
By  a  change  of  variables  r  =  \J2o2v,  and  using  the  definition  of  ?,  it  can  be  shown  that 
(9)  is  equivalent  to 

/OO 

I0(2y/Uq)dv  :=  p(q,  r),  (10) 

2 

where  r  =  The  detection  probability  (10)  is  the  classical  expression  that  can  be  found 
in  [Levanon  1988]  and  [Minkler  and  Minkler  1990]. 


1.3  The  Marcum  Q-Function 


The  Marcum  Q-Function  arises  in  communications  and  radar  signal  processing  problems. 
It  has  a  long  history  in  the  study  of  target  detection  by  pulsed  radars  [Marcum  1950, 
Marcum  1960  and  Marcum  and  Swerling  I960].  As  pointed  out  in  [Simon  1998  and 
Simon  and  Alouini  2003],  it  occurs  in  performance  analysis  related  to  partially  coherent, 
differentiably  coherent  and  noncoherent  communications.  Performance  measures  where 
the  Q-Function  arise  are  error  probabilities  in  transmission  over  fading  channels,  and 
detection  probability  for  code  acquisition  in  direct  sequence  code  division  multiple  access 
systems  [Corazza  and  Ferrari  2002],  A  specific  example  can  be  found  in  [Chiani  1999], 
who  applies  the  Q-Function  to  performance  evaluation  of  noncoherent  and  differentiably 
coherent  detection  of  digital  modulation  over  Nakagami  fading  channels.  In  radar  signal 
processing,  the  Q-Function  appears  as  a  detection  probability.  In  particular,  it  is  related 
to  the  probability  of  detection  on  n  incoherently  integrated  received  signals,  in  a  Gaussian 
clutter  and  noise  environment  [Helstrom  1968,  Nuttall  1975  and  Shnidman  1989]. 

The  standard  Marcum  Q-Function  is  defined  by  the  integral 

/■OO  _  f  x2+a2  \ 

Q(a,j3):=  /  xe  V  2  /  Io(ax)dx.  (11) 

J0 

To  see  the  connection  (11)  has  to  (10),  introduce  the  transformation  a  =  y/2q  and  / 3  = 
y/2 t.  Then  by  introducing  the  transformation  v  =  in  integral  (11),  we  obtain 

Q  (a,/3)  =  Q(V2 s,V2t) 

=  e-<r  /  xe~~Io(\/2^x)dx 

J  y/2 r 

POO 

=  e-<r  /  e~v Io(2y/^v)du.  (12) 
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Hence  it  follows  that  the  single  pulse  probability  of  detection  p(s,  r)  =  Q(v/2?,  \/2r). 

One  of  the  issues  associated  with  (11),  and  in  particular,  the  generalised  Marcum  Q- 
Function  to  be  introduced  in  Section  2,  is  how  to  estimate  it  well.  Also  of  interest  is 
obtaining  upper  and  lower  bounds  on  it.  A  number  of  authors  have  investigated  these 
problems,  including  [Chiani  1999,  Corazza  and  Ferrari  2002,  Simon  1998,  Simon  and 
Alouini  2003  and  Shnidman  1989]. 


1.4  Contributions  of  this  Report 


This  report  outlines  new  mathematical  research  into  the  detection  probability  (10),  and 
its  generalisation  to  n  incoherently  integrated  pulses  in  Gaussian  clutter.  Specifically, 
[Shnidman  1989]  does  not  identify  a  probabilistic  interpretation  to  this  detection  prob¬ 
ability,  which  can  be  used  as  a  mechanism  for  generating  a  Monte  Carlo  estimator  for 
the  generalised  Marcum  Q-Function.  A  new  interpretation  of  the  detection  probability, 
in  Gaussian  clutter  and  noise,  of  a  number  of  incoherently  integrated  received  signals  is 
explored.  In  particular,  it  will  be  shown  that  this  detection  probability  is  equivalent  to  a 
probability  associated  with  two  independent  Poisson  variables.  This  generalises  the  case 
employed  in  [Weinberg  and  Kyprianou  2005]. 

In  addition,  some  new  interesting  properties  of  this  detection  probability  are  explored 
mathematically.  It  is  shown  that  the  generalised  detection  probability  under  investigation 
is  the  solution  to  two  stochastic  Volterra  integral  equation  of  the  second  kind  [Davis  1962 
and  Tricomi  1957].  This  leads  to  useful  representations  that  facilitate  the  construction 
of  bounds  on  the  detection  probability,  and  consequently,  the  Marcum  Q-Function.  We 
will  investigate  the  performance  of  some  new  bounds  on  the  single  pulse  probability  of 
detection  (10). 
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2  Marcum’s  Q-Function  and  a  Poisson 

Connection 

2.1  Generalised  Marcum  Q-Function 


The  generalised  Marcum  Q-Function  of  order  n  €  IN  [Corazza  and  Ferrari  2002,  Nuttall 
1975,  Simon  1998  and  Simon  and  Alouini  2003]  is  defined  by  the  integral 


ro O  _  f  x^+oT  \ 

Qn(a,/3)  =  — — r  /  xne  V  2  J In_1(ax)dx, 
an  1  J/3 


(13) 


where  Ik(x)  is  the  modified  Bessel  function,  of  the  first  kind,  of  order  k  [Bowman  1958], 
which  can  be  defined  by  the  integral  formula 


h(x)  = 


1  r 

2 tt  y_. 


•  —i0\k  —x  sin  6 


(- ie~w)Ke 


d0  = 


l  rTi 
TT  Jo 


cos  6 


cos (k9)d9, 


(14) 


where  i 2  =  — 1.  We  will  be  exploring  this  function  throughout  this  report.  In  particular,  we 
are  interested  in  it  because  of  its  connection  to  radar  detection  probabilities.  As  pointed 
out  in  [Shnidman  1989],  the  detection  probability  of  n  incoherently  integrated  received 
signals  in  a  Gaussian  noise  and  clutter  environment  is  given  by 

(n— 1) 

PniSiT)  =  j  e~^+ni)In^i{2y/unq)dv,  (15) 

when  using  a  quadratic  detector  via  the  Neyman-Pearson  Criterion. 

In  [Weinberg  and  Kyprianou  2005],  (15)  was  investigated  in  the  case  where  n  =  1.  It  was 
shown  that  (15)  for  n  =  1  is  identical  to  the  probability  that  one  Poisson  random  variable, 
representing  a  threshold  variable,  is  less  than  or  equal  to  another  Poisson,  representing 
the  SNR.  We  will  show  that  this  result  can  be  generalised  to  (15). 

It  can  be  shown  [Helstrom  1968,  Nuttall  1975  and  Shnidman  1989]  that  (15)  is  related  to 
(13)  through 

PniSiT)  =  Qn(v/2rK,V/2r).  (16) 


Hence,  any  new  expression  obtained  for  either  one  of  (13)  and  (15)  can  be  applied  to  the 
other,  through  (16). 


2.2  The  Poisson  Connection 


We  now  derive  a  new  probabilistic  expression  for  (15).  The  following  mathematics  is 
included  for  completeness,  and  is  based  upon  the  analysis  in  [Shnidman  1989].  Recall  that 
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the  Bessel  function  of  the  first  kind  of  order  n  £  IN.  denoted  Jn{x),  has  a  series  expansion 
in  terms  of  Gamma  functions 


Jri  (x) 


~  (-l)fc  (|)n+2fc 

fr'o  kw^i  +  k  +  i)’ 


(17) 


[see  Bowman  1958  and  Tsypkin  and  Tsypkin  1988]  and  the  modified  Bessel  function  In(x) 
of  the  first  kind  of  order  n  G  IN  is  related  to  Jn(x)  through  In(x)  =  (—i)n  Jn{ix).  Hence, 
applying  this  to  (17),  and  using  the  fact  that  we  will  only  be  considering  integral  n,  it  is 
not  difficult  to  construct  the  power  series 

oo  xn+2k 

In^  =  ^2n+2kk\{n  +  k)V  (18) 


Hence,  applying  (18),  to  In- i(2y/m7?),  the  detection  probability  (15)  becomes 


Pn(s,r)  = 


(■ n ?)^ 


r°°  n=l  -v  {y/vm)n  1+2fc 

v  2  e  2^  - — - —du 


(n?) 


k= 0 

k  poo 


k\(n  —  1  +  k)\ 


^Qk\(n-l  +  k)\ 


e-vvn~x +kdu. 


(19) 


It  can  be  shown  that 


poo  poo 

e~uun~1+k —  °~r 


roo 

<dv  =  e_T  /  e~v(y  +  T)n~l+kdv. 

Jo 

An  application  of  a  binomial  expansion  to  (u  +  7-)n_1+fc,  transforms  (20)  to 


l  < 


u  +  T)n-1+ke~vdv=  Y, 
1=0 


n~l+k  (n-  \  +  k)\ri 


j! 


Hence,  applying  (21)  to  (19),  we  arrive  at  the  double  series  expansion 

oo/  \k  n—l+k  7 

=  E  jy 

fc=0  1=0 


(20) 


(21) 


(22) 


The  expansion  (22)  is  not  new  to  radar,  but  can  be  found  in  a  number  of  references,  such  as 
[Shnidman  1989].  The  following  interpretation  has  not  been  identified  previously.  Define 
a  pair  of  independent  random  variables  (Ki(n?),  M2(t))  =  (Po(nc),  Po(r)).  We  observe 
that  (22)  can  be  written  in  the  form 


Pn($,r) 


OO 


E 


e  n<’(n?)fc 
k\ 


n—  1+k 

E 


l=o 


e  Trjf 

j! 


^  F[Hi(n?)  =  A;]F[N2(t)  <  n  -  1  +  k] 

k= 0 


F(H2(t)  <n  —  1  +  Hi(n?)). 


(23) 
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This  result  complements  the  corresponding  result  for  n  =  1  in  [Weinberg  and  Kyprianou 
2005], 

From  a  radar  analysis  point  of  view,  (23)  is  an  unexpected  result.  We  can  explain  this 
result  in  a  heuristic  manner  as  follows.  We  can  think  of  the  random  variable  n—  1  +  Mi(n?) 
as  representing  the  signal  to  noise  ratio.  Similarly,  M2 (r)  is  a  random  variable  representing 
the  detection  threshold.  If  we  condition  on  the  threshold  variable,  then  (23)  is  counting 
the  number  of  exceedences  of  a  randomised  threshold  level.  The  larger  the  SNR  c,  the 
more  exceedences  of  the  threshold,  since  ?  is  also  the  mean  value  of  the  associated  Poisson 
distribution.  Also,  Poisson  variables  are  models  for  “rare”  events,  and  detections,  in  some 
circumstances  such  as  low  SNR,  can  be  thought  of  as  such.  Hence  the  form  of  (23)  is 
intuitive  because  it  is  counting  the  number  of  exceedences  of  a  randomised  threshold, 
using  a  model  for  rare  events. 

We  can  also  apply  (23)  to  the  generalised  Marcum  Q-Function.  Note  that,  with  reference 
to  (16),  (23)  implies  that 


Q  n{a,P) 


(24) 


Both  (23)  and  (24)  lead  to  very  simple  Monte  Carlo  estimators  for  these  respective  func¬ 
tions. 
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3  Stochastic  Representation 


This  Section  is  concerned  with  the  derivation  of  a  further  stochastic  representation  of  the 
detection  probability  (15).  In  order  to  derive  this  result,  we  require  an  expression  for  the 
distributional  differences  between  independent  Poisson  random  variables.  Using  Stein’s 
method,  a  new  integral  expression  for  this  is  obtained,  and  is  applied  to  (23). 


3.1  Stein’s  Method  for  Poisson  Approximation 


Stein’s  method  [Barbour  et.  al.  1992,  Chen  1975  and  Stein  1972]  is  a  general  scheme  that 
enables  one  to  obtain  estimates  of  the  rate  of  convergence  of  probability  distributions. 
[Weinberg  2005]  contains  a  detailed  description  of  the  method,  its  history  and  application. 
Our  application  of  it  is  to  construct  an  expression  measuring  the  distributional  difference 
between  two  independent  Poisson  random  variables.  The  application  of  such  an  expression 
is  to  the  construction  of  further  stochastic  representations  of  (23).  In  this  Subsection  we 
present  a  concise  outline  of  Stein’s  method  for  Poisson  approximation.  The  reader  is 
advised  to  consult  [Weinberg  2005]  for  a  more  comprehensive  discussion. 


The  key  idea  behind  Stein’s  method  for  Poisson  approximation  is  to  find  a  solution  to  the 
Stein  equation 

O'  +  1)  -  jg{j)  =  f{j)  -  Po(A )/,  (25) 

where  f(j)  =  A  >  0  and 


n[ieA] 


1  if  j  £  A; 

0  otherwise. 


If  (25)  is  well  defined,  then  for  a  random  variable  W  with  support  the  nonnegative  integers, 


1E[\g{W  +  1)  -  Wg(W)]  =  P(W  €  A)  -  Po(A){H}.  (26) 


Hence,  to  measure  how  well  the  distribution  of  W  is  approximated  by  that  of  a  Poisson 
random  variable,  we  need  to  bound  or  estimate  the  left  hand  side  of  (26). 

The  probabilistic  approach  to  Stein’s  method  [Barbour  et.  al.  1992]  relates  (25)  to  the 
generator  of  an  immigration-death  process,  with  immigration  rate  A,  and  unit  per  capita 
death  rate  [Ross  1983].  The  generator  of  such  a  Markov  process  is 

Ah(j)  =  A [h(j  +  1)  -  h(j)\  -  j[h(j)  -  h(j  -  1)].  (27) 

Here,  we  have  written  g(j  +  1)  =  h(j  +  1)  —  h(j).  The  probabilistic  form  of  the  Stein 
equation  (25)  is 

Ah(j)  =  f(j)  Po(A )/,  (28) 


pOO 

h(j)  =  ~  IE[/(^(t))-Po(A)/]df, 

Jo 


and  has  solution 


(29) 
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where  Zj(t)  is  the  immigration-death  process,  starting  with  j  individuals  [see  Weinberg 
2005].  The  idea,  in  this  context,  is  to  estimate  the  distributional  differences  in  (28)  by 
estimating  the  expectation  of  the  generator  (27),  using  properties  of  the  function  (29)  and 
Markov  process  theory. 

In  the  next  Subsection,  we  apply  this  to  estimate  the  difference  in  distribution  of  two 
independent  Poisson  random  variables. 


3.2  Differences  of  Poisson  Distributions 


We  require  an  expression  for  the  difference  in  distribution  of  W  =  Po (//)  and  V  =  Po(A). 
Hence,  we  simplify  (26).  Note  that 

IE [Wg{W)\  =  f^g(j)e,  ^ 

3=0  J' 

“jA  U  -  1)! 

=  pEfatW  +  l)].  (30) 

Hence,  applying  (30)  to  (26),  it  follows  that 

1E[A g(W  +  1)  -  Wg(W)\  =  (A  -  n)JE [g{W  +  1)].  (31) 

Consequently,  we  need  to  estimate  the  expectation  1E[<7(W  +  1)].  It  is  useful  to  now  use 
the  probabilistic  interpretation  to  Stein’s  method,  in  order  to  simplify  (31).  Hence,  by 
applying  (29),  and  writing  g(j  +  1)  =  h(j  +  1)  -  h(j), 

COO 

g(j  +  1)  =  -  /  nf(Zj+i(t))  -  f(Zj(t))\dt.  (32) 

Jo 

As  in  [Barbour  et.  al.  1992  and  Lindvall  1992],  we  couple  the  two  processes  in  (32)  as 
follows.  The  process  Zj+\{t)  evolves  as  Zj(t)  with  the  addition  of  an  extra  individual,  with 
independent  Exponentially  distributed  lifetime.  Hence  Zj+\(t)  =  Zj(t)  +  Hr^>t],  where  (  is 
an  independent  Exponential  random  variable  with  mean  1.  Then  after  t  >  £,  the  processes 
couple  together.  It  thus  follows  that 

COO 

g(j  +  1)  =  -Ejf  \>t\[f(Zj(t)  +  l)-f{Zj(t))\dt 

coo 

=  -  e_tIE[/(Z?(t)  +  1)  —  f(Zj(t))]dt.  (33) 

Jo 

In  view  of  (31)  and  (33),  what  we  now  require  is  the  distribution  of  the  immigration-death 
process  Zj(t)  when  j  is  a  Poisson  random  variable  with  mean  /r.  The  following  Lemma 
gives  the  required  result: 


Lemma  3.1  For  the  immigration- death  process  Zj(t),  if  W  =  ’Po(fx),  then  Z\y{t)  = 
Po(A  +  (fi  —  A)e_t). 
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Proof: 


Let  Px(s)  =  E[sA]  be  the  probability  generating  function  of  the  random  variable  X.  It 

can  be  shown  that  if  X  =  Po(i'),  then  Px(s )  =  e~u('1~s\  Also,  if  Y  =  Bin (n,p),  then 
Py(s)  =  (1  —p+ps)n.  The  immigration-death  process  Zj(t )  can  be  decomposed  into  a  sum 
of  two  independent  random  variables.  Namely,  it  can  be  shown  that  Zj(t)  =  Xj(t)  +  Zo(t), 

where  Xj(t)  =  Bin(j,  e_<)  and  Zo(t)  =  Po(A(l  —  e_t))  [see  Barbour  et.  al.,  1992],  By 
using  the  independence  of  these  two  variables,  and  by  conditioning  on  W,  it  follows  that 

pZw(t)  =  nsZw(t)]=nnsZw{t)m 

=  E[E[sXw(t)|IT]E[sZo(t)|IP]] 


=  ]E[i3AV(t)(s)]-p2o(t)(s) 


=  E[(l-e-*  +  e-ts)w}e 


A(i-e-*)(i -«) 


=  e-[A+(/i-A)e  ‘K1  -s); 

implying  the  desired  result.  □ 

Since  we  are  interested  in  cumulative  distribution  functions,  we  make  the  choice  of  A  = 
{0, 1, . . . ,  k},  for  some  k  6  IN.  This  is  equivalent  to  the  choice  of  f(j)  =  H[j  <  k\.  Let 
X(t)  =  Zw(t)-  Then  it  follows  that 


roo 

E[g(W  +  l)]  =  -  e_t[P(A(t)  +  l  <  k)  -P(X(t)  <  k)]dt 

J  0 

roo 

=  /  e"fP  (X(t)  =  k)dt.  (34) 

Jo 

Lemma  3.1  implies  the  point  probabilities  of  X(t)  are  given  by 

P(X(t)  =  k)  =  el->+»-».-‘l  lA  +  -  A>e~y .  (35) 

k\ 

Applying  (35)  to  (34),  and  making  the  transformation  v  =  A  +  {/i  —  A)e~f,  and  without 
loss  of  generality,  considering  A  >  fi,  we  obtain 


(36) 

An  application  of  (36)  to  (31)  results  in 

Po(r)/  -  Po(A)/  =  J*  (37) 


10 


DSTO-RR-0304 


Hence  (37)  expresses  the  differences  in  Poisson  cumulative  probability  distributions  as  an 
integral  over  their  respective  means,  of  a  Poisson  distribution  function  whose  mean  is  the 
integration  variable.  The  following  subsection  will  apply  (37)  to  construct  an  alternative 
expression  for  (23). 

Note  that  if  we  introduce  a  random  variable  Q(n,  A)  =  R(/x,  A),  and  an  independent 
Poisson  random  variable  H(^)  =  Po(z/),  then  it  follows  that  (37)  is  equivalent  to 

Po (fj)f  -  Po(A)/  =  (A  -  /i)P[K(0(A,  M))  =  k}.  (38) 

The  random  variable  M(©(A,  //))  in  (38),  a  Poisson  with  a  random  mean,  is  called  a  mixed 
Poisson  distribution. 


3.3  Stochastic  Representation  I 

We  are  now  in  a  position  to  derive  new  expressions  for  the  detection  probability  (23). 
It  is  necessary  to  introduce  a  number  of  random  variables.  Let  Mi(?r?)  =  Po(n?)  and 
H2 (t)  =  Po(r)  be  the  two  independent  Poisson  random  variables  in  the  representation 
(23).  Introduce  random  variables  Hs(?)  =  Po(<j),  H4 (z^)  =  Po(i/)  and  @(r,  ?)  =  R(r  A 
?, r  V  ?).  We  assume  that  both  Hi(n?)  and  ^3(?)  are  independent,  and  when  n  =  1,  it 
is  assumed  they  are  independent  copies  of  each  other  [Lindvall  1992],  Additionally,  it  is 
assumed  that  ^i(n<j),  ^(i')  and  ©(r, f)  are  pairwise  independent. 

Then  by  conditioning  on  Mi(n?), 

pn(?,r)  =  P[H2(t)  <  n  -  1  +  Ki(n?)] 

OO 

=  ^  P[Hi(?r?)  =  P]P[N2(t)  <  n  —  1  +  k] 
k= 0 

=  ^2  P[Ni(nc)  =  k]  ^P[M2(t)  <  n  —  1  +  k\  —  P[Ha(?)  <  b  —  1  +  k]\ 
k= 0  '  ' 

OO 

+  Y2,  =  fc]P[N3(0  <  n  -  1  +  k\.  (39) 

fc=0 

Since,  by  construction,  Mi(n?)  and  ^3(?)  are  independent  and  Poisson  distributed,  we  have 

OO 

^2  P[Ni(n?)  =  ^]P[^3(?)  <  n  -  1  +  k]  =  P[H3(?)  <  n  -  1  +  tti(n?)] 
k= 0 

=  Pn  (?,?),  (40) 

where  we  have  utilised  the  Poisson  expression  (23)  for  the  detection  probability  (15).  Also, 
with  an  application  of  the  Poisson  difference  equation  (37),  we  have 

PDM'r)  <  n  —  1  +  k]  —  P[^3(?)  <  n  —  1  +  k] 
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fS  e~f vn-l+k 


k\ 


dv 


=  J  P[N4(z/)  =  n  —  1  +  k]du.  (41) 

Consequently,  by  applying  (41)  and  using  the  fact  that  Hi  (n?)  and  K4(i/)  are  independent, 
f]P[Hi(n?)  =  k\  (p[H2(t)  <  n-  1  +  k\  -P[H3(?)  <  n-  1  +  k]\ 

k= 0  '  ' 

=  j  P[K4(z/)  =i  n  —  1  +  ^i(nq)]du.  (42) 

Finally,  combining  (40)  and  (42),  (39)  becomes 

Pn{^,r)  =  p„(s,?)  +  J  P[H4(^)  =  n  -  1  +  Hi(n?)]dzA  (43) 

This  gives  a  new  expression  for  the  detection  probability  (15).  It  shows  that  the  latter  is 
equal  to  the  same  detection  probability,  where  the  threshold  is  equal  to  the  SNR,  plus  or 
minus  a  discrepancy  factor  that  depends  on  the  difference  between  r  and  ?. 


It  is  interesting  to  note  that  the  integrand  in  (43)  is  closely  related  to  the  detection 
probability  pn(?,  u)  =  P[H4(z/)  <  n  —  1  +  Ki(n<t)].  In  fact,  it  is  obvious  that  P[H4(z/)  = 
n  —  1  +  Hi(n?)]  =  pn(s,v)  —  P[K4(z/)  <  n  —  1  +  Hi(n?)].  Hence,  the  integrand  is  a 
nonlinear  function  of  pn{<>,  v)-  This  suggests  that  the  represention  (43)  is  a  nonlinear 
integral  equation  of  pn(?,  v).  as  a  function  of  its  second  independent  variable,  namely  v. 

As  pointed  out  in  [Davis  1962  and  Tricomi  1957],  a  function  f(x),  defined  on  a  suitable 
domain,  is  the  solution  to  a  nonlinear  Volterra  integral  equation  if  it  satisfies 

f(x)  =  f(x o)+  [  F[t,f(t)\dt,  (44) 

j  X  0 

where  F  is  nonlinear.  In  the  case  of  Volterra  integral  equations  of  the  first  and  second 
kind,  it  is  assumed  that  the  function  F  is  a  product  of  the  function  /  and  a  kernel  K.  We 
illustrate  how  (43)  is  of  the  form  (44). 


Write  the  integrand  in  (43)  as 


P[H4(i/)  =  n  —  1  +  Ki(n?)]  =  pn(^,u)Kn(q,v), 
where  the  kernel  Kn(s,  u)  is 


P[N4(i/)  =  n  —  1  +  Ki(n?)] 
~~  P[H4(z/)  <  n-  1  +  Ni(n?)]‘ 


(45) 


(46) 


Consider  ?  (and  n)  hxed,  and  choose  f(x)  =  pn(s,x),  F[x,f(x)]  =  —f(x)Kn($,x)  and 
xq  =  C-  Then  (43)  is  exactly  of  the  form  (44),  as  a  function  of  the  threshold,  namely 
x  =  T. 
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It  is  worth  noting  that  a  function  /  satisfying  the  integral  equation  (44)  can  be  esti¬ 
mated  through  a  Picard  iteration  scheme  [Giles  1987].  One  considers  the  functional  series 
{/m,m  €  IN}  with  fo('x)  =  f(x q)  and  fm(x)  defined  recursively  through 

fm{x)  =  f(x0)  +  f  F[t,  fm-i(t)\dt.  (47) 

Jx  0 

There  are  existence  and  uniqueness  theorems,  which  give  conditions  on  F  to  guaran¬ 
tee  a  reasonable  approximation  is  obtained  [Giles  1987  and  Tricomi  1957].  In  principle, 
one  could  apply  the  scheme  (47)  to  (43)  to  obtain  estimates  of  the  detection  probability 
pn(s,  t).  In  addition,  it  may  be  possible  to  use  properties  of  the  solution  to  Volterra  inte¬ 
gral  equations  to  derive  new  bounds  on  the  Marcum  Q-function.  The  validity  and  merit 
of  these  approaches  will  be  investigated  in  future  research. 

Observe  that,  by  taking  partial  derivatives,  with  respect  to  r,  it  is  not  difficult  to  see  that 

J^PnObr)  =  -P[N4(t)  =  n  -  1  +  Ki(n?)],  (48) 

implying  that  pn(s,  r)  is  a  decreasing  function  of  the  threshold  r,  for  fixed  ?.  Note  that, 
in  view  of  the  preceding  comments,  the  partial  derivative  (48)  is  a  nonlinear  function  of 
PniSi  T)- 


It  is  still  possible  to  write  (43)  in  a  more  compact  stochastic  form.  Using  the  definition  of 
0(t,  c),  it  is  not  difficult  to  see  that 


Pn(?,r)  =  Pn(?,?)  +  (?  -  r)  P[H4(0(r,  ?))  =  n  -  1  +  Ki(n?)].  (49) 

Both  (43)  and  (49)  are  new  expressions  for  the  detection  probability  (15),  and  both  can  be 
applied  to  the  generalised  Marcum  Q-Function  (13),  by  an  application  of  (16),  to  derive 
analogous  expressions. 


3.4  Stochastic  Representation  II 

It  is  possible  to  derive  a  representation  analogous  to  (43),  by  instead  conditioning  on 
^2(7"),  instead  of  conditioning  on  Mi(n<j).  Introduce  a  random  variable  M5 (jit)  =  Po(nr), 
which  we  assume  to  be  independent  of  M2  (r).  Throughout  we  use  the  same  definitions  of 
appropriate  random  variables,  as  defined  in  Subsection  3.3.  Then,  using  a  similar  argument 
as  in  the  derivation  of  (43), 

Pnis,r)  =  P[K2(t)  <  n  -  1  +  Mi(n?)] 


OO 

=  ^2  P[^2 (t" )  =  k]  (P[Hi(n?)  >  k  —  n  +  1]  —  P^s^rr)  >  k  —  n  +  1]) 
k= 0 

OO 

+  ^2  P[^2 (7" )  =  &]P[N5(nr)  >  k  —  n  +  1] 
k= 0 
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=  ^  P[K2(t)  =  k]  (P[^5(nr)  <  k  —  n\  —  PfN^ns)  <  k  —  n}) 

k=0 

+P[H2(r)  <  H5  (nr)  +  n  —  1] 

f-n<; 

=  +  /  P[H4(z/)  =  H2(t)  -  n\du 

Jnr 

=  Pu{t,t)  +  J  nP[tt4(?rz/)  =  H2(t)  -  n]diA  (50) 

Note  that,  analogously  to  (48),  by  taking  partial  derivatives  of  (50)  with  respect  to 

^Pnis.r)  =  nP[M4(n?)  =  K2(t)  -n].  (51) 

Consequently,  (51)  shows  that  pn(s,  r)  is  an  increasing  function  of  the  SNR  parameter  ?, 
for  fixed  r. 

We  can  also  construct  an  analogue  of  (49).  Using  the  same  definitions  as  in  the  latter 
expression, 

Pn(s,T)  =  Pn{r,T )  +  n(q  -  r)P[N4(n0(r,  ?))  =  H2(r)  -n].  (52) 

The  representations  (43)  and  (50)  indicate  that  the  detection  probability  (15),  and  also 
the  Marcum  Q-Function  (13),  may  be  estimated  using  properties  of  Volterra  integral 
equations.  This  idea  is  partially  explored  in  the  next  Section. 
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4  Bounds  on  the  Pulse  Detection  Probability 


In  this  final  Section  we  consider  bounding  the  detection  probability  (15)  in  the  simple  case 
of  n  =  1.  We  will  refer  to  this  case  as  the  pulse  detection  probability.  Bounding  this  case, 
and  equivalently,  obtaining  bounds  on  the  standard  Marcum  Q-Function  (11),  has  been 
the  subject  of  much  interest  in  recent  years  [Chiani  1999,  Corazza  and  Ferrari  2002,  Simon 
1998  and  Simon  and  Alouini  (2000  and  2003)].  We  investigate  whether  bounds  derived 
from  the  new  stochastic  representations  (43)  and  (50)  can  improve  known  bounds  on  this 
pulse  detection  probability. 


It  is  important  to  note  that  a  closed  form  result  exists  for  the  pulse  detection  probability, 
in  the  case  where  ?  =  r.  Specifically,  it  can  be  shown  [Schwartz,  Bennett  and  Stein  1966] 


that 


1  +  e-2M0(2?) 


(53) 


This  result  will  be  an  integral  component  of  the  new  bounds.  A  probabilistic  proof  of  (53) 
can  be  found  in  Appendix  A.  Before  deriving  some  new  bounds  based  on  (43)  and  (50), 
we  examine  some  well-known  bounds. 


4.1  Lower  and  Upper  Bounds 

To  begin,  we  examine  known  lower  and  upper  bounds  of  the  detection  probability  (15)  in 
the  case  of  n  =  1.  Such  bounds  have  been  examined  for  the  Marcum  Q-Function  (13), 
for  n  =  1,  in  [Chiani  1999,  Corazza  and  Ferrari  2002,  Simon  1998  and  Simon  and  Alouini 
(2000  and  2003)].  In  the  following  we  will  express  their  bounds  in  terms  of  the  single  pulse 
probability  of  detection.  Tight  lower  bounds  give  estimates  of  the  least  possible  values  of 
the  detection  probability,  and  thus  have  been  used  as  a  performance  measure  in  [Chiani 
1999].  Upper  bounds  indicate  the  maximum  possible  detection  probability,  and  are  also 
of  interest  in  performance  analysis.  Unless  otherwise  stated,  it  will  be  understood  that  all 
lower  bounds  are  taken,  in  practice,  as  the  maximum  of  the  expressed  bound  and  zero,  in 
order  to  avoid  meaningless  lower  bounds.  Similarly,  all  upper  bounds  are  the  minimum  of 
the  expressed  bound  and  unity. 


4.1.1  Case  1 :  q  >  r 


An  upper  bound  on  the  pulse  detection  probability,  or  equivalently,  the  standard  Marcum 
Q-function,  for  the  case  where  ^  >  r  have  only  recently  appeared  in  [Corazza  and  Ferrari 
2002],  These  authors  derive  the  upper  bound 


p(ct) 


<  !  Jo(2V^)  \  -c  c-UU2WC2?i2 

eW  \ 


+  [Erfc(— y^)  -  Erfc(Vr  -  y7?)]  j 


(54) 
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where  Erfc  is  the  complementary  error  function 

2  r°°  •> 

Erfc (2:)  :=  —=  /  dt.  (55) 

V7T  Jz 

As  [Corazza  and  Ferrari  2002]  emphasize,  the  upper  bound  (54)  is  the  only  known  non¬ 
trivial  upper  bound  on  the  pulse  detection  probability  for  the  case  under  consideration. 

Lower  bounds  for  ?  >  r  are  more  prevalent  in  the  literature.  A  lower  bound  analogous  to 
(54)  is  also  introduced  in  [Corazza  and  Ferrari  2002],  which  is 

p(q,T)  >  1  -  e-|[2^2]  |e-#  -  e-W'fc-tf 


Wf 


Erfc  I -4)-  Erfc  (ALA 

y/2J  l  V2 


(56) 


where  7?  =  log(Jo6V^))  _ 


A  Chernoff-like  lower  bound  is  derived  in  [Simon  and  Alouini  2000],  given  by 


p(S,t)  >  1  - 


1 


=-[C?-C?]2  _e-[^+xA]2 


[Chiani  1999]  also  proposes  a  Chernoff-like  lower  bound,  which  is 

p(s,t)  >  e~^+T]I0(2^r). 

Finally,  [Simon  1998]  suggests  the  bound 

1 


p(s,t)  >  1 


,-k-T]2 


(57) 


(58) 


(59) 


We  now  include  some  comments  on  the  performance  and  accuracy  of  these  bounds.  It  is 
reported  in  [Corazza  and  Ferrari  2002]  that  the  upper  bound  (54)  performs  well,  when 
compared  to  the  exact  pulse  detection  probability,  except  as  r  increases.  [Chiani  1999] 
states  that  the  lower  bound  (57)  is  very  tight  and  better  than  (58).  [Corazza  and  Ferrari 
2002]  also  report  that  their  lower  bound  (56)  is  also  very  tight,  and  only  matched  by  the 
bound  (57).  The  bound  (59)  is  reported  to  be  quite  unreliable  by  [Corazza  and  Ferrari 
2002], 


4.1.2  Case  2:  q  <  t 


In  this  case,  there  are  many  proposed  upper  and  lower  bounds  in  the  literature.  [Corazza 
and  Ferrari  2002]  derive  the  two  sided  bound 


/o(2yAr) 

g2v/?r 


v/vrrErfc(v/T:  -  y/q)  <  p(q,  r) 

Io(2v/?r) 


-  e2C?F  ^ e 


je  [r  ?l2  +  0r?Er ic(y/r  -  y7?)}  •  (60) 
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[Simon  and  Alouini  2000]  present  the  bound 

g-h+d2  <  p(<j, r)  <  e-!'1"-*'!2.  (61) 

[Chiani  1999]  derives  the  two  sided  bound 

g-h+r]  j0(2v/?r)  <  p(q,  t)  <  e_^+T^I0(2v/?r)  +  ^Erfc(V^  -  V?).  (62) 

The  final  bound  we  consider  is  that  due  to  [Simon  1998],  which  is 

/l  re -[T+?]2  <  P<S,t)  <  ^  e-^12.  (63) 

V^+  y/T- 


Again,  we  provide  some  comments  on  the  accuracy  and  performance  of  these  bounds. 
[Corazza  and  Ferrari  2002]  compare  their  upper  bound  with  the  three  others  introduced 
here,  and  found  that  their  upper  bound  in  (60)  has  the  best  performance.  Their  lower 
bound  in  (60)  also  performs  well,  and  is  only  matched  by  that  in  (62).  [Chiani  1999] 
reports  that  the  bounds  in  (62)  are  tighter  than  those  in  (61). 


4.2  Some  New  Lower  and  Upper  Bounds 


In  this  section  we  derive  some  new  upper  and  lower  bounds  on  the  detection  probability 
(15)  in  the  case  where  n  =  1.  These  bounds  are  derived  using  the  representations  (43) 
and  (50).  We  consider  the  two  cases  ?  >  r  and  q  <  r  separately.  Observe  that  in  view  of 
(43)  and  (50),  we  have 


=  p(q,q)+ j  P[Ni (i/)  =  bf2(?)]^,  (64) 

and 

p(q,r)  =  p{r,r)  +  J  IP[^i(^)  —  N2(t)  —  l]du,  (65) 

where  Hi(-)  and  ^(O  are  independent  Poisson  random  variables,  with  appropriate  means. 
Throughout  we  will  utilise  the  notation  that  a  V  b  is  the  maximum  of  a  and  b ,  while  a  A  6 
is  the  minimum  of  a  and  b. 


4.2.1  Case  1:  q  >  r 

Note  that  in  this  case  both  integrals  in  (64)  and  (65)  are  nonnegative,  and  so  an  appropriate 
lower  bound  is  given  by 

pis,  t )  >  p{q ,  q)  V  p(r,  r),  (66) 

where  we  can  apply  (53)  to  evaluate  this  exactly. 
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In  order  to  derive  an  upper  bound,  recall  that  (51)  implies  p{q,r)  is  an  increasing  function 
in  q.  Hence, 

p(c,t)  <  + j  F[^2(?)  <  {v)\dv 

=  p(s,q)  + j  p(y,s)dv 

<  [l  +  ?-r]p(?,?).  (67) 

By  applying  a  similar  argument  to  (65),  and  using  the  fact  that  in  view  of  (48),  p(q,r)  is 
a  decreasing  function  of  r,  it  can  be  shown  that 

p(s,t)  <  [1  +  q -t]p(t,t).  (68) 

Hence,  by  combining  (67)  and  (68),  we  arrive  at  the  bound 

Pis,  r)  <  [1  +  ?  -  r]  (p(q,  q)  A  p(r,  r)) .  (69) 

Finally,  by  combining  the  bounds  (66)  and  (69),  we  arrive  at  the  two  sided  bound 

(p(<b  0  V  p(r,  t))  <  p(q,  t)  <  [1  +  ?  —  t]  (p(q,  q)  A  p(r,  r)) .  (70) 

On  inspection  of  the  upper  bound  (69)  it  is  clear  that  if  the  difference  q  —  r  is  much  larger 
than  minimum  of  the  two  probabilities,  then  the  bound  will  be  likely  to  exceed  unity. 

4.2.2  Case  2:  q  <  r 

In  this  case,  we  note  that  the  integrals  in  (64)  and  (65)  are  nonpositive,  and  so  it  follows 
that  we  have  the  upper  bound 

p(q,r)  <  p(q,q)  A  p(t,t).  (71) 

To  derive  a  lower  bound,  note  that 

IP[Hi(z/)  =  tt2(?)]  <  p( v,  q)  <  p{q,  ?),  (72) 

and  so  applying  (72)  to  (64),  we  arrive  at 

p{q,T)>[l  +  q-T]p(q,q).  (73) 

Using  a  similar  argument  applied  to  (65),  it  can  be  shown  that 

p(q,r)  >[l  +  q -t]p(t,t).  (74) 

Hence,  by  combining  (73)  and  (74),  we  arrive  at 

p(?,  r)  >  (1  +  ?  -  r)(p(s,  q)  V  p(r,  t)).  (75) 

Thus,  by  combining  (71)  and  (75),  we  arrive  at  the  two  sided  bound 

(1  +  ?  -  r)(p(?,  q)  V  p(r,  r))  <  p(?,  r)  <  (p(?,  ?)  A  p(r,  r))  (76) 

Note  that  the  lower  bound  (75)  may  also  experience  the  same  problems  as  the  upper  bound 
(69),  resulting  in  a  trivial  lower  bound  in  (76). 
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4.2.3  Some  Comments  on  the  Representation  in  Equation  (64) 

Consider  again  the  representation  of  (64),  which  we  write  in  the  form 

p(S,r)  =  +  J  P[Ni(i/)  -K2(?)  =  0}du,  (77) 

where  ^i(-)  and  H2O)  are  independent  Poisson  random  variables.  Define  the  random 
variable 

Z(i/,?)  =  N1(i/)-N2(?).  (78) 

The  difference  of  independent  Poisson  random  variables  has  been  studied  extensively  in 
the  literature  [Irwin  1937,  Karlis  and  Ntzoufras  2003  and  Skellam  1946],  and  is  known  as  a 
Skellam  distribution.  This  distribution  has  found  application  in  the  analysis  of  sports  data 
[see  Karlis  and  Ntzoufras  2003  and  references  contained  therein].  Thus,  the  probability 
in  the  integral  component  of  (77)  is  the  probability  that  a  Skellam  distribution  takes  the 
value  zero.  Hence  it  may  be  possible  to  find  some  good  approximations,  and  bounds,  for 
W[Z(y,  <j)  =  0]  in  the  literature.  A  preliminary  review  has  found  Gaussian  approximations 
for  this,  but  they  were  considered  not  to  be  useful  for  the  applications  considered  in  this 
report.  Part  of  the  problem  is  that  the  random  variable  in  (78)  has  mean  v  —  and 

variance  zz+?,  so  that  a  Gaussian  approximation  applied  to  P[— 1  <  Z(u,  ?)  <  1]  will  result 
in  intractible  integrals  in  (77).  What  is  needed  is  an  easily  integrated  approximation. 

In  further  research  in  this  area,  the  author  plans  to  examine  the  literature  more  extensively, 
for  work  on  the  Skellam  distribution,  in  order  to  improve  the  bounds  of  Subsection  4.2. 

In  the  next  Section  we  perform  some  numerical  comparisons  to  investigate  whether  the 
new  bounds  (70)  and  (76)  provide  any  improvement  on  the  bounds  in  Subsection  4.1. 


4.3  Numerical  Comparison  of  Bounds 

We  are  now  in  a  position  to  compare  the  bounds  of  Subsection  4.2  to  those  in  Subsection 
4.1.  Extensive  numerical  experiments  showed  that  the  new  bounds  of  Subsection  4.2  can 
provide  some  improvements  on  the  corresponding  bounds  in  Subsection  4.1,  but  are  not 
globally  better.  We  consider  a  selection  of  cases  to  show  the  strengths  and  weaknesses  of 
these  new  bounds. 

Each  figure  includes  a  comparison  of  the  bounds  to  an  almost  exact  result.  The  latter  has 
been  obtained  via  recursive  adaptive  Simpson  quadrature,  with  a  tolerance  of  10~6. 

Figures  1  and  2  are  plots  of  upper  bounds  for  the  case  where  c  >  t.  In  this  case,  we 
compare  the  upper  bound  (54),  due  to  [Corazza  and  Ferrari  2002]  to  the  new  upper  bound 
in  (70).  Note  that  the  vertical  scale  is  in  natural  logarithms.  Each  Figure  contains  two 
subplots.  The  first  subplot  shows  the  bounds  as  a  function  of  r,  for  a  fixed  ?,  while  the 
second  subplot  shows  it  as  a  function  of  ?,  for  a  fixed  r.  The  legend  used  in  the  Figure 
refers  to  the  bound  of  [Corazza  and  Ferrari  2002]  as  ‘C-F’,  while  Wei’  refers  to  the  bound 
in  (70).  The  probability  estimated  from  adaptive  Simpson  quadrature  is  referred  to  as 
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‘Exact’.  As  can  be  observed,  the  new  upper  bound  provides  no  improvement  on  that  due 
to  [Corazza  and  Ferrari  2002],  Other  numerical  experiments  provided  no  improvement  on 
this  result.  This  is  most  likely  due  to  the  fact  that  more  work  needs  to  be  done  to  obtain 
a  satisfactory  estimate  of  the  integral  component  of  (77). 


Comparison  of  Upper  Bounds  for  g  >  x,  for  fixed  g 


Comparison  of  Upper  Bounds  for  g  >  x,  for  fixed  x 


Figure  1:  Comparison  of  upper  bounds  for  the  case  where  g  >  t,  showing  the  bounds  (54) 
and  (70)  as  functions  of  r  and  g.  ‘C-F’  refers  to  the  upper  bound  due  to  Corazza  and 
Ferrari,  namely  (54),  while  ‘Wei’  is  the  upper  bound  in  (70).  The  vertical  scale  is  in 
logarithms. 

We  now  consider  upper  bounds  for  the  case  where  g  <  r.  In  this  situation,  we  have  four 
upper  bounds  to  compare  to  the  new  upper  bound  in  (76).  Specifically,  we  are  interested 
in  how  the  upper  bounds  in  (60)-(63)  compare  to  (76).  Two  sets  of  (g,r)  parameters  are 
considered.  Figures  3  and  4  are  for  one  case,  while  Figures  5  and  6  are  a  different  case. 
Bounds  in  the  Figures  are  abbreviated  as  ‘C-F’,  for  that  due  to  [Corazza  and  Ferrari  2002], 
‘S-A’,  for  [Simon  and  Alouini  2000],  ‘Sim’  for  [Simon  1998],  ‘Chi’  for  [Chiani  1999]  and 
‘Wei’  for  that  in  (76).  As  before,  the  vertical  scale  is  in  logarithms.  The  first  subplot  of 
Figure  3  shows  that  the  new  bound  (76)  performs  better  than  all  others,  with  the  exception 
of  the  upper  bound  (62),  due  to  [Chiani  1999].  The  second  subplot  shows  a  case  where 
the  new  bound  has  superior  performance.  Figure  4  is  the  same  as  Figure  3,  except  the 
two  worst  performing  bounds  have  been  removed. 

Figures  5  and  6  show  the  upper  bounds  (60)-(63)  compared  to  (76),  in  a  slightly  different 
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Comparison  of  Upper  Bounds  for  q  >  x,  for  fixed  q 
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Figure  2:  Similar  plot  to  that  of  Figure  1,  except  on  a  different  range  of  q  and  r  values. 


scenario.  In  this  example,  the  upper  bound  due  to  [Chiani  1999]  has  the  best  performance, 
while  the  upper  bound  in  (76)  tends  to  be  very  inaccurate.  Figure  6  is  the  same  as  Figure 
5,  with  the  removal  of  the  two  worst  performing  bounds. 

We  now  investigate  lower  bounds,  beginning  with  the  case  where  q  >  r.  Figures  7  and 
8  are  for  one  particular  selection  of  (<j,  r)  parameters.  The  lower  bounds  investigated  are 
those  in  (56)- (59),  and  the  new  lower  bound  in  (70).  The  bounds  are  described  in  the 
legend  using  the  same  abbreviations  as  before.  Although  difficult  to  see,  the  exact  values 
lie  almost  on  the  horizontal  axis.  Figures  7  and  8  shows  that  the  new  bound  (70)  performs 
well,  and  is  only  outperformed  by  the  lower  bound  (56)  due  to  [Corazza  and  Ferrari  2002], 

The  final  scenario  we  consider  is  lower  bounds  in  the  case  where  q  <  r.  Figure  9  shows  a 
plot  of  the  lower  bound  components  of  (60)- (63),  as  well  as  the  new  lower  bound  in  (76). 
The  same  naming  conventions  are  employed  in  the  legend,  and  the  vertical  scale  is  also  in 
logarithms.  The  first  subplot  of  Figure  9  is  an  example  where  the  new  lower  bound  in  (76) 
has  superior  performance.  The  second  subplot  in  Figure  9  shows  the  same  effect,  except 
the  lower  bound  in  (62)  is  slightly  better  for  small  values  of  r.  Figure  10  is  the  same  as 
Figure  9,  except  the  two  best  performing  bounds  are  compared. 
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Comparison  of  Upper  Bounds  for  g  <  x,  for  fixed  g 


Comparison  of  Upper  Bounds  for  g  <  x,  for  fixed  x 
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Figure  3:  A  comparison  of  upper  bounds  in  the  case  where  <g  <  r.  The  new  upper  bound 
in  (76),  labelled  ‘Wei’,  is  outperformed  by  the  bound  (62)  due  to  Chiani,  labelled  as  ‘Chi’. 
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Comparison  of  Upper  Bounds  for  q  <  x,  for  fixed  q 


X 

Comparison  of  Upper  Bounds  for  q  <  x,  for  fixed  x 


Figure  4 ■'  The  same  plot  as  for  Figure  3,  with  the  removal  of  the  two  worst  performing 
bounds. 
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Comparison  of  Upper  Bounds  for  g  <  x,  for  fixed  g 


Comparison  of  Upper  Bounds  for  g  <  x,  for  fixed  x 


Figure  5:  Comparison  of  upper  bounds  in  the  case  where  g  <  r,  but  over  a  different  range 
of  g  and  r  parameters.  In  both  subplots,  it  is  clear  that  the  bounds  ‘Sim  due  to  Simon 
(63)  and  ‘S-A’,  due  to  Simon  and  Alouini  (61)  are  unsatisfactory.  Figure  6  shows  the 
same  plot  with  these  two  upper  bounds  removed. 


24 


DSTO-RR-0304 


Comparison  of  Upper  Bounds  for  g  <  x,  for  fixed  g 


Comparison  of  Upper  Bounds  for  g  <  x,  for  fixed  x 


Figure  6:  This  is  the  same  as  Figure  5,  except  the  two  worst  performing  bounds  have  been 
removed.  It  is  now  clear  that  the  bound  labelled  ‘C-F’,  due  to  Corazza  and  Ferrari  (60)  is 
also  unsatisfactory.  The  bound  due  to  Chiani  (62)  performs  the  best,  while  the  new  upper 
bound  in  (76),  labelled  as  ‘Wei’,  is  the  second  best. 
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Comparison  of  Lower  Bounds  for  q  >  x,  for  fixed  q 


Comparison  of  Lower  Bounds  for  q  >  x,  for  fixed  x 


Figure  7:  Plots  of  lower  bounds  in  the  case  where  (,  >  t.  The  same  labelling  convention 
has  been  employed  as  in  previous  figures.  The  exact  value,  labelled  as  ‘Exact’,  is  almost 
on  the  horizontal  axis  in  each  subplot.  It  is  clear  that  the  new  bound  (70)  is  performing 
very  well.  See  Figure  8  for  clarification  of  this. 


26 


DSTO-RR-0304 


Comparison  of  Lower  Bounds  for  q  >  x,  for  fixed  q 


Comparison  of  Lower  Bounds  for  q  >  x,  for  fixed  x 


Figure  8:  The  same  plot  as  for  Figure  7,  except  with  the  removal  of  the  bound  (59)  due  to 
Simon.  As  for  Figure  7,  the  exact  value  is  almost  on  the  horizontal  axis.  The  best  bound 
is  that  due  to  Corazza  and  Ferrari  (56),  while  the  new  lower  bound  in  (70),  namely  ‘Wei’, 
is  second  best. 
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Figure  9:  Lower  bounds  for  the  case  where  g  <  r.  This  is  an  example  where  the  new  lower 
bound  in  (76)  has  very  good  performance.  The  only  bound  that  is  comparable  to  the  new 
lower  bound  is  that  due  to  Chiani  (62).  The  latter  is  better  in  the  second  subplot  when  g 
is  very  small. 
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Comparison  of  Lower  Bounds  for  q  <  x,  for  fixed  q 


Comparison  of  Lower  Bounds  for  q  <  x,  for  fixed  x 
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Figure  1 0:  The  same  plot  as  in  Figure  9,  with  the  inclusion  of  the  Chiani  lower  bound  in 
(62),  and  the  new  lower  bound  in  (76),  and  the  exclusion  of  the  other  lower  bounds.  The 
first  subplot  shows  the  improvement  more  clearly,  while  the  second  subplot  shows  there  is 
a  point  where  the  better  bound  switches  from  the  Chiani  one  to  the  new  lower  bound. 
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5  Conclusions  and  Future  Directions 

This  report  introduced  a  number  of  new  results.  A  Poisson  association  with  the  Mar¬ 
cum  Q-Function,  and  associated  detection  probabilities,  was  derived,  extending  the  case 
considered  in  [Weinberg  and  Kyprianou  2005].  Stein’s  method  was  used  to  produce  a 
new  expression  for  the  distributional  differences  of  a  pair  of  Poisson  random  variables. 
This  permitted  the  construction  of  two  Volterra  integral  equations,  whose  solution  was 
the  detection  probability  under  consideration.  These  new  representations  were  used  to 
construct  lower  and  upper  bounds  on  the  single  pulse  probability  of  detection.  Numerical 
comparisons  to  other  bounds  in  the  literature  showed  there  are  cases  where  these  new 
bounds  provided  improvements. 

It  is  believed  that  better  global  lower  and  upper  bounds  could  be  achieved  through  further 
analysis  of  these  Volterra  integral  equations.  Specifically,  there  may  be  properties  known 
about  such  equations  that  will  lead  to  new  ways  of  estimating  the  Marcum  Q-Function, 
and  associated  detection  probabilities.  As  pointed  out  in  the  report,  it  may  be  possible 
to  estimate  the  detection  probability  (15)  using  a  Picard  iteration  scheme  applied  to  (43). 
Also,  using  estimates  of  the  Skellam  distribution’s  zero  probability,  applied  to  (77),  may 
also  produce  improved  bounds.  This  will  be  the  focus  of  subsequent  research  in  this  area. 
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Appendix  A:  Proof  of  Equation  (53) 


Suppose  X  and  Y  are  independent  and  identically  distributed  Poisson  random  variables 
with  mean  A.  Then,  with  reference  to  the  expansion  (18), 

°o  /  -A \k\2  °o  \2k 

v{x=Y)=h\—)  =e"“Sw=e"M/o(2A)'  (A-1) 

Since  X  and  Y  are  statistically  identical,  it  follows  that 

P(X  <Y)  =  P(X  >  Y).  (A. 2) 


Thus,  (A. 2)  implies  that 

P(X  +  Y)  =  2JP(X  <  Y).  (A. 3) 

Hence,  with  an  application  of  (A.l)  and  (A. 3),  it  follows  that 


PpT<F)  =  P(X  =  Y)  +  P(X  <  Y) 


=  e~2XI0(  2A)  +  \[l~  e~2XI0(  2\j 

=  |  [l  +  e-2AT0(2A)]  .  (A. 4) 

The  desired  result  now  follows  by  an  application  of  (A.4)  to  (23)  in  the  case  where  n  =  1. 
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